home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 21 / AACD 21.iso / AACD / Utilities / Ghostscript / src / gscie.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-01-01  |  38.2 KB  |  1,278 lines
  1. /* Copyright (C) 1992, 1995, 1996, 1997, 1998, 1999 Aladdin Enterprises.  All rights reserved.
  2.   
  3.   This file is part of AFPL Ghostscript.
  4.   
  5.   AFPL Ghostscript is distributed with NO WARRANTY OF ANY KIND.  No author or
  6.   distributor accepts any responsibility for the consequences of using it, or
  7.   for whether it serves any particular purpose or works at all, unless he or
  8.   she says so in writing.  Refer to the Aladdin Free Public License (the
  9.   "License") for full details.
  10.   
  11.   Every copy of AFPL Ghostscript must include a copy of the License, normally
  12.   in a plain ASCII text file named PUBLIC.  The License grants you the right
  13.   to copy, modify and redistribute AFPL Ghostscript, but only under certain
  14.   conditions described in the License.  Among other things, the License
  15.   requires that the copyright notice and this notice be preserved on all
  16.   copies.
  17. */
  18.  
  19. /*$Id: gscie.c,v 1.3 2000/09/19 19:00:26 lpd Exp $ */
  20. /* CIE color rendering cache management */
  21. #include "math_.h"
  22. #include "memory_.h"
  23. #include "gx.h"
  24. #include "gserrors.h"
  25. #include "gsstruct.h"
  26. #include "gsmatrix.h"        /* for gscolor2.h */
  27. #include "gxcspace.h"        /* for gxcie.c */
  28. #include "gscolor2.h"        /* for gs_set/currentcolorrendering */
  29. #include "gxarith.h"
  30. #include "gxcie.h"
  31. #include "gxdevice.h"        /* for gxcmap.h */
  32. #include "gxcmap.h"
  33. #include "gzstate.h"
  34.  
  35. /* Forward references */
  36. private int cie_joint_caches_init(P3(gx_cie_joint_caches *,
  37.                      const gs_cie_common *,
  38.                      gs_cie_render *));
  39. private void cie_joint_caches_complete(P4(gx_cie_joint_caches *,
  40.                       const gs_cie_common *,
  41.                       const gs_cie_abc *,
  42.                       const gs_cie_render *));
  43. private void cie_cache_restrict(P2(cie_cache_floats *, const gs_range *));
  44. private void cie_mult3(P3(const gs_vector3 *, const gs_matrix3 *,
  45.               gs_vector3 *));
  46. private void cie_matrix_mult3(P3(const gs_matrix3 *, const gs_matrix3 *,
  47.                  gs_matrix3 *));
  48. private void cie_invert3(P2(const gs_matrix3 *, gs_matrix3 *));
  49. private void cie_matrix_init(P1(gs_matrix3 *));
  50.  
  51. /* Allocator structure types */
  52. private_st_joint_caches();
  53.  
  54. #define RESTRICTED_INDEX(v, n, itemp)\
  55.   ((uint)(itemp = (int)(v)) >= (n) ?\
  56.    (itemp < 0 ? 0 : (n) - 1) : itemp)
  57.  
  58. /* Define the template for loading a cache. */
  59. /* If we had parameterized types, or a more flexible type system, */
  60. /* this could be done with a single procedure. */
  61. #define CIE_LOAD_CACHE_BODY(pcache, domains, rprocs, dprocs, pcie, cname)\
  62.   BEGIN\
  63.     int j;\
  64. \
  65.     for (j = 0; j < countof(pcache); j++) {\
  66.         cie_cache_floats *pcf = &(pcache)[j].floats;\
  67.         int i;\
  68.         gs_for_loop_params lp;\
  69. \
  70.         gs_cie_cache_init(&pcf->params, &lp, &(domains)[j], cname);\
  71.         for (i = 0; i < gx_cie_cache_size; lp.init += lp.step, i++) {\
  72.         pcf->values[i] = (*(rprocs)->procs[j])(lp.init, pcie);\
  73.         if_debug5('C', "[C]%s[%d,%d] = %g => %g\n",\
  74.               cname, j, i, lp.init, pcf->values[i]);\
  75.         }\
  76.         pcf->params.is_identity =\
  77.         (rprocs)->procs[j] == (dprocs).procs[j];\
  78.     }\
  79.   END
  80.  
  81. /*
  82.  * Determine whether a function is a linear transformation of the form
  83.  * f(x) = scale * x + origin.
  84.  */
  85. private bool
  86. cache_is_linear(cie_linear_params_t *params, const cie_cache_floats *pcf)
  87. {
  88.     double origin = pcf->values[0];
  89.     double diff = pcf->values[countof(pcf->values) - 1] - origin;
  90.     double scale = diff / (countof(pcf->values) - 1);
  91.     int i;
  92.     double test = origin + scale;
  93.  
  94.     for (i = 1; i < countof(pcf->values) - 1; ++i, test += scale)
  95.     if (fabs(pcf->values[i] - test) >= 0.5 / countof(pcf->values))
  96.         return (params->is_linear = false);
  97.     params->origin = origin - pcf->params.base;
  98.     params->scale = diff * pcf->params.factor / (countof(pcf->values) - 1);
  99.     return (params->is_linear = true);
  100. }
  101.  
  102. private void
  103. cache_set_linear(cie_cache_floats *pcf)
  104. {
  105.     if (pcf->params.is_identity) {
  106.         if_debug1('c', "[c]linear(0x%lx) = true (is_identity)\n",
  107.               (ulong)pcf);
  108.         pcf->params.linear.is_linear = true;
  109.         pcf->params.linear.origin = 0;
  110.         pcf->params.linear.scale = 1;
  111.     } else if (cache_is_linear(&pcf->params.linear, pcf)) {
  112.         if (pcf->params.linear.origin == 0 &&
  113.         fabs(pcf->params.linear.scale - 1) < 0.00001)
  114.         pcf->params.is_identity = true;
  115.         if_debug4('c',
  116.               "[c]linear(0x%lx) = true, origin = %g, scale = %g%s\n",
  117.               (ulong)pcf, pcf->params.linear.origin,
  118.               pcf->params.linear.scale,
  119.               (pcf->params.is_identity ? " (=> is_identity)" : ""));
  120.     }
  121. #ifdef DEBUG
  122.     else
  123.         if_debug1('c', "[c]linear(0x%lx) = false\n", (ulong)pcf);
  124. #endif
  125. }
  126. private void
  127. cache3_set_linear(gx_cie_vector_cache *caches /*[3]*/)
  128. {
  129.     cache_set_linear(&caches[0].floats);
  130.     cache_set_linear(&caches[1].floats);
  131.     cache_set_linear(&caches[2].floats);
  132. }
  133.  
  134. #ifdef DEBUG
  135. private void
  136. if_debug_vector3(const char *str, const gs_vector3 *vec)
  137. {
  138.     if_debug4('c', "%s[%g %g %g]\n", str, vec->u, vec->v, vec->w);
  139. }
  140. private void
  141. if_debug_matrix3(const char *str, const gs_matrix3 *mat)
  142. {
  143.     if_debug10('c', "%s [%g %g %g] [%g %g %g] [%g %g %g]\n", str,
  144.            mat->cu.u, mat->cu.v, mat->cu.w,
  145.            mat->cv.u, mat->cv.v, mat->cv.w,
  146.            mat->cw.u, mat->cw.v, mat->cw.w);
  147. }
  148. #else
  149. #  define if_debug_vector3(str, vec) DO_NOTHING
  150. #  define if_debug_matrix3(str, mat) DO_NOTHING
  151. #endif
  152.  
  153. /* ------ Default values for CIE dictionary elements ------ */
  154.  
  155. /* Default transformation procedures. */
  156.  
  157. private float
  158. a_identity(floatp in, const gs_cie_a * pcie)
  159. {
  160.     return in;
  161. }
  162. private float
  163. a_from_cache(floatp in, const gs_cie_a * pcie)
  164. {
  165.     return gs_cie_cached_value(in, &pcie->caches.DecodeA.floats);
  166. }
  167.  
  168. private float
  169. abc_identity(floatp in, const gs_cie_abc * pcie)
  170. {
  171.     return in;
  172. }
  173. private float
  174. abc_from_cache_0(floatp in, const gs_cie_abc * pcie)
  175. {
  176.     return gs_cie_cached_value(in, &pcie->caches.DecodeABC[0].floats);
  177. }
  178. private float
  179. abc_from_cache_1(floatp in, const gs_cie_abc * pcie)
  180. {
  181.     return gs_cie_cached_value(in, &pcie->caches.DecodeABC[1].floats);
  182. }
  183. private float
  184. abc_from_cache_2(floatp in, const gs_cie_abc * pcie)
  185. {
  186.     return gs_cie_cached_value(in, &pcie->caches.DecodeABC[2].floats);
  187. }
  188.  
  189. private float
  190. def_identity(floatp in, const gs_cie_def * pcie)
  191. {
  192.     return in;
  193. }
  194. private float
  195. def_from_cache_0(floatp in, const gs_cie_def * pcie)
  196. {
  197.     return gs_cie_cached_value(in, &pcie->caches_def.DecodeDEF[0].floats);
  198. }
  199. private float
  200. def_from_cache_1(floatp in, const gs_cie_def * pcie)
  201. {
  202.     return gs_cie_cached_value(in, &pcie->caches_def.DecodeDEF[1].floats);
  203. }
  204. private float
  205. def_from_cache_2(floatp in, const gs_cie_def * pcie)
  206. {
  207.     return gs_cie_cached_value(in, &pcie->caches_def.DecodeDEF[2].floats);
  208. }
  209.  
  210. private float
  211. defg_identity(floatp in, const gs_cie_defg * pcie)
  212. {
  213.     return in;
  214. }
  215. private float
  216. defg_from_cache_0(floatp in, const gs_cie_defg * pcie)
  217. {
  218.     return gs_cie_cached_value(in, &pcie->caches_defg.DecodeDEFG[0].floats);
  219. }
  220. private float
  221. defg_from_cache_1(floatp in, const gs_cie_defg * pcie)
  222. {
  223.     return gs_cie_cached_value(in, &pcie->caches_defg.DecodeDEFG[1].floats);
  224. }
  225. private float
  226. defg_from_cache_2(floatp in, const gs_cie_defg * pcie)
  227. {
  228.     return gs_cie_cached_value(in, &pcie->caches_defg.DecodeDEFG[2].floats);
  229. }
  230. private float
  231. defg_from_cache_3(floatp in, const gs_cie_defg * pcie)
  232. {
  233.     return gs_cie_cached_value(in, &pcie->caches_defg.DecodeDEFG[3].floats);
  234. }
  235.  
  236. private float
  237. common_identity(floatp in, const gs_cie_common * pcie)
  238. {
  239.     return in;
  240. }
  241. private float
  242. lmn_from_cache_0(floatp in, const gs_cie_common * pcie)
  243. {
  244.     return gs_cie_cached_value(in, &pcie->caches.DecodeLMN[0].floats);
  245. }
  246. private float
  247. lmn_from_cache_1(floatp in, const gs_cie_common * pcie)
  248. {
  249.     return gs_cie_cached_value(in, &pcie->caches.DecodeLMN[1].floats);
  250. }
  251. private float
  252. lmn_from_cache_2(floatp in, const gs_cie_common * pcie)
  253. {
  254.     return gs_cie_cached_value(in, &pcie->caches.DecodeLMN[2].floats);
  255. }
  256.  
  257. /* Transformation procedures for accessing an already-loaded cache. */
  258.  
  259. float
  260. gs_cie_cached_value(floatp in, const cie_cache_floats *pcache)
  261. {
  262.     /*
  263.      * We need to get the same results when we sample an already-loaded
  264.      * cache, so we need to round the index just a tiny bit.
  265.      */
  266.     int index =
  267.     (int)((in - pcache->params.base) * pcache->params.factor + 0.0001);
  268.  
  269.     CIE_CLAMP_INDEX(index);
  270.     return pcache->values[index];
  271. }
  272.  
  273. /* Default vectors and matrices. */
  274.  
  275. const gs_range3 Range3_default = {
  276.     { {0, 1}, {0, 1}, {0, 1} }
  277. };
  278. const gs_range4 Range4_default = {
  279.     { {0, 1}, {0, 1}, {0, 1}, {0, 1} }
  280. };
  281. const gs_cie_defg_proc4 DecodeDEFG_default = {
  282.     {defg_identity, defg_identity, defg_identity, defg_identity}
  283. };
  284. const gs_cie_defg_proc4 DecodeDEFG_from_cache = {
  285.     {defg_from_cache_0, defg_from_cache_1, defg_from_cache_2, defg_from_cache_3}
  286. };
  287. const gs_cie_def_proc3 DecodeDEF_default = {
  288.     {def_identity, def_identity, def_identity}
  289. };
  290. const gs_cie_def_proc3 DecodeDEF_from_cache = {
  291.     {def_from_cache_0, def_from_cache_1, def_from_cache_2}
  292. };
  293. const gs_cie_abc_proc3 DecodeABC_default = {
  294.     {abc_identity, abc_identity, abc_identity}
  295. };
  296. const gs_cie_abc_proc3 DecodeABC_from_cache = {
  297.     {abc_from_cache_0, abc_from_cache_1, abc_from_cache_2}
  298. };
  299. const gs_cie_common_proc3 DecodeLMN_default = {
  300.     {common_identity, common_identity, common_identity}
  301. };
  302. const gs_cie_common_proc3 DecodeLMN_from_cache = {
  303.     {lmn_from_cache_0, lmn_from_cache_1, lmn_from_cache_2}
  304. };
  305. const gs_matrix3 Matrix3_default = {
  306.     {1, 0, 0},
  307.     {0, 1, 0},
  308.     {0, 0, 1},
  309.     1 /*true */
  310. };
  311. const gs_range RangeA_default = {0, 1};
  312. const gs_cie_a_proc DecodeA_default = a_identity;
  313. const gs_cie_a_proc DecodeA_from_cache = a_from_cache;
  314. const gs_vector3 MatrixA_default = {1, 1, 1};
  315. const gs_vector3 BlackPoint_default = {0, 0, 0};
  316.  
  317. /* Initialize a CIE color. */
  318. /* This only happens on setcolorspace. */
  319. void
  320. gx_init_CIE(gs_client_color * pcc, const gs_color_space * pcs)
  321. {
  322.     gx_init_paint_4(pcc, pcs);
  323.     /* (0...) may not be within the range of allowable values. */
  324.     (*pcs->type->restrict_color)(pcc, pcs);
  325. }
  326.  
  327. /* Restrict CIE colors. */
  328.  
  329. inline private void
  330. cie_restrict(float *pv, const gs_range *range)
  331. {
  332.     if (*pv <= range->rmin)
  333.     *pv = range->rmin;
  334.     else if (*pv >= range->rmax)
  335.     *pv = range->rmax;
  336. }
  337.  
  338. void
  339. gx_restrict_CIEDEFG(gs_client_color * pcc, const gs_color_space * pcs)
  340. {
  341.     const gs_cie_defg *pcie = pcs->params.defg;
  342.  
  343.     cie_restrict(&pcc->paint.values[0], &pcie->RangeDEFG.ranges[0]);
  344.     cie_restrict(&pcc->paint.values[1], &pcie->RangeDEFG.ranges[1]);
  345.     cie_restrict(&pcc->paint.values[2], &pcie->RangeDEFG.ranges[2]);
  346.     cie_restrict(&pcc->paint.values[3], &pcie->RangeDEFG.ranges[3]);
  347. }
  348. void
  349. gx_restrict_CIEDEF(gs_client_color * pcc, const gs_color_space * pcs)
  350. {
  351.     const gs_cie_def *pcie = pcs->params.def;
  352.  
  353.     cie_restrict(&pcc->paint.values[0], &pcie->RangeDEF.ranges[0]);
  354.     cie_restrict(&pcc->paint.values[1], &pcie->RangeDEF.ranges[1]);
  355.     cie_restrict(&pcc->paint.values[2], &pcie->RangeDEF.ranges[2]);
  356. }
  357. void
  358. gx_restrict_CIEABC(gs_client_color * pcc, const gs_color_space * pcs)
  359. {
  360.     const gs_cie_abc *pcie = pcs->params.abc;
  361.  
  362.     cie_restrict(&pcc->paint.values[0], &pcie->RangeABC.ranges[0]);
  363.     cie_restrict(&pcc->paint.values[1], &pcie->RangeABC.ranges[1]);
  364.     cie_restrict(&pcc->paint.values[2], &pcie->RangeABC.ranges[2]);
  365. }
  366. void
  367. gx_restrict_CIEA(gs_client_color * pcc, const gs_color_space * pcs)
  368. {
  369.     const gs_cie_a *pcie = pcs->params.a;
  370.  
  371.     cie_restrict(&pcc->paint.values[0], &pcie->RangeA);
  372. }
  373.  
  374. /* ================ Table setup ================ */
  375.  
  376. /* ------ Install a CIE color space ------ */
  377.  
  378. private void cie_load_common_cache(P2(gs_cie_common *, gs_state *));
  379. private void cie_cache_mult(P3(gx_cie_vector_cache *, const gs_vector3 *,
  380.                    const cie_cache_floats *));
  381. private bool cie_cache_mult3(P2(gx_cie_vector_cache *,
  382.                 const gs_matrix3 *));
  383.  
  384. private int
  385. gx_install_cie_abc(gs_cie_abc *pcie, gs_state * pgs)
  386. {
  387.     if_debug_matrix3("[c]CIE MatrixABC =", &pcie->MatrixABC);
  388.     cie_matrix_init(&pcie->MatrixABC);
  389.     CIE_LOAD_CACHE_BODY(pcie->caches.DecodeABC, pcie->RangeABC.ranges,
  390.             &pcie->DecodeABC, DecodeABC_default, pcie,
  391.             "DecodeABC");
  392.     cie_load_common_cache(&pcie->common, pgs);
  393.     gs_cie_abc_complete(pcie);
  394.     return gs_cie_cs_complete(pgs, true);
  395. }
  396.  
  397. int
  398. gx_install_CIEDEFG(const gs_color_space * pcs, gs_state * pgs)
  399. {
  400.     gs_cie_defg *pcie = pcs->params.defg;
  401.  
  402.     CIE_LOAD_CACHE_BODY(pcie->caches_defg.DecodeDEFG, pcie->RangeDEFG.ranges,
  403.             &pcie->DecodeDEFG, DecodeDEFG_default, pcie,
  404.             "DecodeDEFG");
  405.     return gx_install_cie_abc((gs_cie_abc *)pcie, pgs);
  406. }
  407.  
  408. int
  409. gx_install_CIEDEF(const gs_color_space * pcs, gs_state * pgs)
  410. {
  411.     gs_cie_def *pcie = pcs->params.def;
  412.  
  413.     CIE_LOAD_CACHE_BODY(pcie->caches_def.DecodeDEF, pcie->RangeDEF.ranges,
  414.             &pcie->DecodeDEF, DecodeDEF_default, pcie,
  415.             "DecodeDEF");
  416.     return gx_install_cie_abc((gs_cie_abc *)pcie, pgs);
  417. }
  418.  
  419. int
  420. gx_install_CIEABC(const gs_color_space * pcs, gs_state * pgs)
  421. {
  422.     return gx_install_cie_abc(pcs->params.abc, pgs);
  423. }
  424.  
  425. int
  426. gx_install_CIEA(const gs_color_space * pcs, gs_state * pgs)
  427. {
  428.     gs_cie_a *pcie = pcs->params.a;
  429.     int i;
  430.     gs_for_loop_params lp;
  431.     float in;
  432.  
  433.     gs_cie_cache_init(&pcie->caches.DecodeA.floats.params, &lp,
  434.               &pcie->RangeA, "DecodeA");
  435.     for (i = 0, in = lp.init; i < gx_cie_cache_size; in += lp.step, i++) {
  436.     pcie->caches.DecodeA.floats.values[i] = (*pcie->DecodeA)(in, pcie);
  437.     if_debug3('C', "[C]DecodeA[%d] = %g => %g\n",
  438.           i, in, pcie->caches.DecodeA.floats.values[i]);
  439.     }
  440.     cie_load_common_cache(&pcie->common, pgs);
  441.     gs_cie_a_complete(pcie);
  442.     return gs_cie_cs_complete(pgs, true);
  443. }
  444.  
  445. /* Load the common caches when installing the color space. */
  446. private void
  447. cie_load_common_cache(gs_cie_common * pcie, gs_state * pgs)
  448. {
  449.     if_debug_matrix3("[c]CIE MatrixLMN =", &pcie->MatrixLMN);
  450.     cie_matrix_init(&pcie->MatrixLMN);
  451.     CIE_LOAD_CACHE_BODY(pcie->caches.DecodeLMN, pcie->RangeLMN.ranges,
  452.             &pcie->DecodeLMN, DecodeLMN_default, pcie,
  453.             "DecodeLMN");
  454. }
  455.  
  456. /* Complete loading the common caches. */
  457. private void
  458. cie_common_complete(gs_cie_common *pcie)
  459. {
  460.     int i;
  461.  
  462.     for (i = 0; i < 3; ++i)
  463.     cache_set_linear(&pcie->caches.DecodeLMN[i].floats);
  464. }
  465.  
  466. /*
  467.  * Restrict the DecodeDEF[G] cache according to RangeHIJ[K], and scale to
  468.  * the dimensions of Table.
  469.  */
  470. private void
  471. gs_cie_defx_scale(float *values, const gs_range *range, int dim)
  472. {
  473.     double scale = (dim - 1.0) / (range->rmax - range->rmin);
  474.     int i;
  475.  
  476.     for (i = 0; i < gx_cie_cache_size; ++i) {
  477.     float value = values[i];
  478.  
  479.     values[i] =
  480.         (value <= range->rmin ? 0 :
  481.          value >= range->rmax ? dim - 1 :
  482.          (value - range->rmin) * scale);
  483.     }
  484. }
  485.  
  486. /* Complete loading a CIEBasedDEFG color space. */
  487. /* This routine is NOT idempotent. */
  488. void
  489. gs_cie_defg_complete(gs_cie_defg * pcie)
  490. {
  491.     int j;
  492.  
  493.     for (j = 0; j < 4; ++j)
  494.     gs_cie_defx_scale(pcie->caches_defg.DecodeDEFG[j].floats.values,
  495.               &pcie->RangeHIJK.ranges[j], pcie->Table.dims[j]);
  496.     gs_cie_abc_complete((gs_cie_abc *)pcie);
  497. }
  498.  
  499. /* Complete loading a CIEBasedDEF color space. */
  500. /* This routine is NOT idempotent. */
  501. void
  502. gs_cie_def_complete(gs_cie_def * pcie)
  503. {
  504.     int j;
  505.  
  506.     for (j = 0; j < 3; ++j)
  507.     gs_cie_defx_scale(pcie->caches_def.DecodeDEF[j].floats.values,
  508.               &pcie->RangeHIJ.ranges[j], pcie->Table.dims[j]);
  509.     gs_cie_abc_complete((gs_cie_abc *)pcie);
  510. }
  511.  
  512. /* Complete loading a CIEBasedABC color space. */
  513. /* This routine is idempotent. */
  514. void
  515. gs_cie_abc_complete(gs_cie_abc * pcie)
  516. {
  517.     cache3_set_linear(pcie->caches.DecodeABC);
  518.     pcie->caches.skipABC =
  519.     cie_cache_mult3(pcie->caches.DecodeABC, &pcie->MatrixABC);
  520.     cie_common_complete((gs_cie_common *)pcie);
  521. }
  522.  
  523. /* Complete loading a CIEBasedA color space. */
  524. /* This routine is idempotent. */
  525. void
  526. gs_cie_a_complete(gs_cie_a * pcie)
  527. {
  528.     cie_cache_mult(&pcie->caches.DecodeA, &pcie->MatrixA,
  529.            &pcie->caches.DecodeA.floats);
  530.     cache_set_linear(&pcie->caches.DecodeA.floats);
  531.     cie_common_complete((gs_cie_common *)pcie);
  532. }
  533.  
  534. /* Convert a scalar cache to a vector cache by multiplying */
  535. /* the scalar values by a vector. */
  536. /* This procedure is idempotent. */
  537. private void
  538. cie_cache_mult(gx_cie_vector_cache * pcache, const gs_vector3 * pvec,
  539.            const cie_cache_floats * pcf)
  540. {
  541.     int i;
  542.  
  543.     pcache->vecs.params.base = float2cie_cached(pcf->params.base);
  544.     pcache->vecs.params.factor = float2cie_cached(pcf->params.factor);
  545.     pcache->vecs.params.limit =
  546.     float2cie_cached((gx_cie_cache_size - 1) / pcf->params.factor +
  547.              pcf->params.base);
  548.     for (i = 0; i < gx_cie_cache_size; ++i) {
  549.     float f = pcf->values[i];
  550.  
  551.     pcache->vecs.values[i].u = float2cie_cached(f * pvec->u);
  552.     pcache->vecs.values[i].v = float2cie_cached(f * pvec->v);
  553.     pcache->vecs.values[i].w = float2cie_cached(f * pvec->w);
  554.     }
  555. }
  556.  
  557. /* Convert 3 scalar caches to vector caches by multiplying by a matrix. */
  558. /* Return true iff the resulting cache is an identity transformation. */
  559. private bool
  560. cie_cache_mult3(gx_cie_vector_cache * pc /*[3] */ , const gs_matrix3 * pmat)
  561. {
  562.     cie_cache_mult(pc, &pmat->cu, &pc->floats);
  563.     cie_cache_mult(pc + 1, &pmat->cv, &pc[1].floats);
  564.     cie_cache_mult(pc + 2, &pmat->cw, &pc[2].floats);
  565.     return pmat->is_identity & pc[0].floats.params.is_identity &
  566.     pc[1].floats.params.is_identity & pc[2].floats.params.is_identity;
  567. }
  568.  
  569. /* ------ Install a rendering dictionary ------ */
  570.  
  571. /* setcolorrendering */
  572. int
  573. gs_setcolorrendering(gs_state * pgs, gs_cie_render * pcrd)
  574. {
  575.     int code = gs_cie_render_complete(pcrd);
  576.     const gs_cie_render *pcrd_old = pgs->cie_render;
  577.     bool joint_ok;
  578.  
  579.     if (code < 0)
  580.     return code;
  581.     if (pcrd_old != 0 && pcrd->id == pcrd_old->id)
  582.     return 0;        /* detect needless reselecting */
  583.     joint_ok =
  584.     pcrd_old != 0 &&
  585. #define CRD_SAME(elt) !memcmp(&pcrd->elt, &pcrd_old->elt, sizeof(pcrd->elt))
  586.     CRD_SAME(points.WhitePoint) && CRD_SAME(points.BlackPoint) &&
  587.     CRD_SAME(MatrixPQR) && CRD_SAME(RangePQR) &&
  588.     CRD_SAME(TransformPQR);
  589. #undef CRD_SAME
  590.     rc_assign(pgs->cie_render, pcrd, "gs_setcolorrendering");
  591.     /* Initialize the joint caches if needed. */
  592.     if (!joint_ok)
  593.     code = gs_cie_cs_complete(pgs, true);
  594.     gx_unset_dev_color(pgs);
  595.     return code;
  596. }
  597.  
  598. /* currentcolorrendering */
  599. const gs_cie_render *
  600. gs_currentcolorrendering(const gs_state * pgs)
  601. {
  602.     return pgs->cie_render;
  603. }
  604.  
  605. /* Unshare (allocating if necessary) the joint caches. */
  606. gx_cie_joint_caches *
  607. gx_currentciecaches(gs_state * pgs)
  608. {
  609.     gx_cie_joint_caches *pjc = pgs->cie_joint_caches;
  610.  
  611.     rc_unshare_struct(pgs->cie_joint_caches, gx_cie_joint_caches,
  612.               &st_joint_caches, pgs->memory,
  613.               return 0, "gx_currentciecaches");
  614.     if (pgs->cie_joint_caches != pjc) {
  615.     pjc = pgs->cie_joint_caches;
  616.     pjc->cspace_id = pjc->render_id = gs_no_id;
  617.     pjc->id_status = pjc->status = CIE_JC_STATUS_BUILT;
  618.     }
  619.     return pjc;
  620. }
  621.  
  622. /* Compute the parameters for loading a cache, setting base and factor. */
  623. /* This procedure is idempotent. */
  624. void
  625. gs_cie_cache_init(cie_cache_params * pcache, gs_for_loop_params * pflp,
  626.           const gs_range * domain, client_name_t cname)
  627. {
  628.     /* We need to map the values in the range
  629.      * [domain->rmin..domain->rmax].  However, if rmin < 0 < rmax and
  630.      * the function is non-linear, this can lead to anomalies at zero,
  631.      * which is the default value for CIE colors.  The "correct" way to
  632.      * approach this is to run the mapping functions on demand, but we
  633.      * don't want to deal with the complexities of the callbacks this
  634.      * would involve (especially in the middle of rendering images);
  635.      * instead, we adjust the range so that zero maps precisely to a
  636.      * cache slot.  Define:
  637.      *      a = domain->rmin;
  638.      *      b = domain->rmax;
  639.      *      R = b - a;
  640.      *      N = gx_cie_cache_size - 1;
  641.      *      f(v) = N(v-a)/R;
  642.      *      x = f(0).
  643.      * If x is not an integer, we can either increase b or
  644.      * decrease a to make it one.  In the former case, compute:
  645.      *      Kb = floor(x); R'b = N(0-a)/Kb; b' = a + R'b.
  646.      * In the latter case, compute:
  647.      *      Ka = ceiling(x-N); R'a = N(0-b)/Ka; a' = b - R'a.
  648.      * We choose whichever method stretches the range the least,
  649.      * i.e., the one whose R' value (R'a or R'b) is smaller.
  650.      */
  651.     double a = domain->rmin, b = domain->rmax;
  652.     double R = b - a;
  653.     double delta;
  654. #define N (gx_cie_cache_size - 1)
  655.  
  656.     /* Adjust the range if necessary. */
  657.     if (a < 0 && b >= 0) {
  658.     double x = -N * a / R;    /* must be > 0 */
  659.     double Kb = floor(x);    /* must be >= 0 */
  660.     double Ka = ceil(x) - N;    /* must be <= 0 */
  661.  
  662.     if (Kb == 0 || (Ka != 0 && -b / Ka < -a / Kb))    /* use R'a */
  663.         R = -N * b / Ka, a = b - R;
  664.     else            /* use R'b */
  665.         R = -N * a / Kb, b = a + R;
  666.     }
  667.     delta = R / N;
  668. #ifdef CIE_CACHE_INTERPOLATE
  669.     pcache->base = a;        /* no rounding */
  670. #else
  671.     pcache->base = a - delta / 2;    /* so lookup will round */
  672. #endif
  673.     pcache->factor = (delta == 0 ? 0 : N / R);
  674.     if_debug4('c', "[c]cache %s 0x%lx base=%g, factor=%g\n",
  675.           (const char *)cname, (ulong) pcache,
  676.           pcache->base, pcache->factor);
  677.     pflp->init = a;
  678.     pflp->step = delta;
  679.     pflp->limit = b + delta / 2;
  680. #undef N
  681. }
  682.  
  683. /* ------ Complete a rendering structure ------ */
  684.  
  685. /*
  686.  * Compute the derived values in a CRD that don't involve the cached
  687.  * procedure values.  This procedure is idempotent.
  688.  */
  689. private void cie_transform_range3(P3(const gs_range3 *, const gs_matrix3 *,
  690.                      gs_range3 *));
  691. int
  692. gs_cie_render_init(gs_cie_render * pcrd)
  693. {
  694.     gs_matrix3 PQR_inverse;
  695.  
  696.     if (pcrd->status >= CIE_RENDER_STATUS_INITED)
  697.     return 0;        /* init already done */
  698.     if_debug_matrix3("[c]CRD MatrixLMN =", &pcrd->MatrixLMN);
  699.     cie_matrix_init(&pcrd->MatrixLMN);
  700.     if_debug_matrix3("[c]CRD MatrixABC =", &pcrd->MatrixABC);
  701.     cie_matrix_init(&pcrd->MatrixABC);
  702.     if_debug_matrix3("[c]CRD MatrixPQR =", &pcrd->MatrixPQR);
  703.     cie_matrix_init(&pcrd->MatrixPQR);
  704.     cie_invert3(&pcrd->MatrixPQR, &PQR_inverse);
  705.     cie_matrix_mult3(&pcrd->MatrixLMN, &PQR_inverse,
  706.              &pcrd->MatrixPQR_inverse_LMN);
  707.     cie_transform_range3(&pcrd->RangePQR, &pcrd->MatrixPQR_inverse_LMN,
  708.              &pcrd->DomainLMN);
  709.     cie_transform_range3(&pcrd->RangeLMN, &pcrd->MatrixABC,
  710.              &pcrd->DomainABC);
  711.     cie_mult3(&pcrd->points.WhitePoint, &pcrd->MatrixPQR, &pcrd->wdpqr);
  712.     cie_mult3(&pcrd->points.BlackPoint, &pcrd->MatrixPQR, &pcrd->bdpqr);
  713.     pcrd->status = CIE_RENDER_STATUS_INITED;
  714.     return 0;
  715. }
  716.  
  717. /*
  718.  * Sample the EncodeLMN, EncodeABC, and RenderTableT CRD procedures, and
  719.  * load the caches.  This procedure is idempotent.
  720.  */
  721. int
  722. gs_cie_render_sample(gs_cie_render * pcrd)
  723. {
  724.     int code;
  725.  
  726.     if (pcrd->status >= CIE_RENDER_STATUS_SAMPLED)
  727.     return 0;        /* sampling already done */
  728.     code = gs_cie_render_init(pcrd);
  729.     if (code < 0)
  730.     return code;
  731.     CIE_LOAD_CACHE_BODY(pcrd->caches.EncodeLMN, pcrd->DomainLMN.ranges,
  732.             &pcrd->EncodeLMN, Encode_default, pcrd, "EncodeLMN");
  733.     cache3_set_linear(pcrd->caches.EncodeLMN);
  734.     CIE_LOAD_CACHE_BODY(pcrd->caches.EncodeABC, pcrd->DomainABC.ranges,
  735.             &pcrd->EncodeABC, Encode_default, pcrd, "EncodeABC");
  736.     if (pcrd->RenderTable.lookup.table != 0) {
  737.     int i, j, m = pcrd->RenderTable.lookup.m;
  738.     gs_for_loop_params flp;
  739.     bool is_identity = true;
  740.  
  741.     for (j = 0; j < m; j++) {
  742.         gs_cie_cache_init(&pcrd->caches.RenderTableT[j].fracs.params,
  743.                   &flp, &Range3_default.ranges[0],
  744.                   "RenderTableT");
  745.         is_identity &= pcrd->RenderTable.T.procs[j] ==
  746.         RenderTableT_default.procs[j];
  747.     }
  748.     pcrd->caches.RenderTableT_is_identity = is_identity;
  749.     /*
  750.      * Unfortunately, we defined the first argument of the RenderTable
  751.      * T procedures as being a byte, limiting the number of distinct
  752.      * cache entries to 256 rather than gx_cie_cache_size.
  753.      * We confine this decision to this loop, rather than propagating
  754.      * it to the procedures that use the cached data, so that we can
  755.      * change it more easily at some future time.
  756.      */
  757.     for (i = 0; i < gx_cie_cache_size; i++) {
  758. #if gx_cie_log2_cache_size >= 8
  759.         byte value = i >> (gx_cie_log2_cache_size - 8);
  760. #else
  761.         byte value = (i << (8 - gx_cie_log2_cache_size)) +
  762.         (i >> (gx_cie_log2_cache_size * 2 - 8));
  763. #endif
  764.         for (j = 0; j < m; j++) {
  765.         pcrd->caches.RenderTableT[j].fracs.values[i] =
  766.             (*pcrd->RenderTable.T.procs[j])(value, pcrd);
  767.         if_debug3('C', "[C]RenderTableT[%d,%d] = %g\n",
  768.               i, j,
  769.               frac2float(pcrd->caches.RenderTableT[j].fracs.values[i]));
  770.         }
  771.     }
  772.     }
  773.     pcrd->status = CIE_RENDER_STATUS_SAMPLED;
  774.     return 0;
  775. }
  776.  
  777. /* Transform a set of ranges. */
  778. private void
  779. cie_transform_range(const gs_range3 * in, floatp mu, floatp mv, floatp mw,
  780.             gs_range * out)
  781. {
  782.     float umin = mu * in->ranges[0].rmin, umax = mu * in->ranges[0].rmax;
  783.     float vmin = mv * in->ranges[1].rmin, vmax = mv * in->ranges[1].rmax;
  784.     float wmin = mw * in->ranges[2].rmin, wmax = mw * in->ranges[2].rmax;
  785.     float temp;
  786.  
  787.     if (umin > umax)
  788.     temp = umin, umin = umax, umax = temp;
  789.     if (vmin > vmax)
  790.     temp = vmin, vmin = vmax, vmax = temp;
  791.     if (wmin > wmax)
  792.     temp = wmin, wmin = wmax, wmax = temp;
  793.     out->rmin = umin + vmin + wmin;
  794.     out->rmax = umax + vmax + wmax;
  795. }
  796. private void
  797. cie_transform_range3(const gs_range3 * in, const gs_matrix3 * mat,
  798.              gs_range3 * out)
  799. {
  800.     cie_transform_range(in, mat->cu.u, mat->cv.u, mat->cw.u,
  801.             &out->ranges[0]);
  802.     cie_transform_range(in, mat->cu.v, mat->cv.v, mat->cw.v,
  803.             &out->ranges[1]);
  804.     cie_transform_range(in, mat->cu.w, mat->cv.w, mat->cw.w,
  805.             &out->ranges[2]);
  806. }
  807.  
  808. /*
  809.  * Finish preparing a CRD for installation, by restricting and/or
  810.  * transforming the cached procedure values.
  811.  * This procedure is idempotent.
  812.  */
  813. int
  814. gs_cie_render_complete(gs_cie_render * pcrd)
  815. {
  816.     int code;
  817.  
  818.     if (pcrd->status >= CIE_RENDER_STATUS_COMPLETED)
  819.     return 0;        /* completion already done */
  820.     code = gs_cie_render_sample(pcrd);
  821.     if (code < 0)
  822.     return code;
  823.     /*
  824.      * Since range restriction happens immediately after
  825.      * the cache lookup, we can save a step by restricting
  826.      * the values in the cache entries.
  827.      *
  828.      * If there is no lookup table, we want the final ABC values
  829.      * to be fracs; if there is a table, we want them to be
  830.      * appropriately scaled ints.
  831.      */
  832.     pcrd->MatrixABCEncode = pcrd->MatrixABC;
  833.     {
  834.     int c;
  835.     double f;
  836.  
  837.     for (c = 0; c < 3; c++) {
  838.         gx_cie_float_fixed_cache *pcache = &pcrd->caches.EncodeABC[c];
  839.  
  840.         cie_cache_restrict(&pcrd->caches.EncodeLMN[c].floats,
  841.                    &pcrd->RangeLMN.ranges[c]);
  842.         cie_cache_restrict(&pcrd->caches.EncodeABC[c].floats,
  843.                    &pcrd->RangeABC.ranges[c]);
  844.         if (pcrd->RenderTable.lookup.table == 0) {
  845.         cie_cache_restrict(&pcache->floats,
  846.                    &Range3_default.ranges[0]);
  847.         gs_cie_cache_to_fracs(&pcache->floats, &pcache->fixeds.fracs);
  848.         pcache->fixeds.fracs.params.is_identity = false;
  849.         } else {
  850.         int i;
  851.         int n = pcrd->RenderTable.lookup.dims[c];
  852.  
  853. #ifdef CIE_RENDER_TABLE_INTERPOLATE
  854. #  define SCALED_INDEX(f, n, itemp)\
  855.      RESTRICTED_INDEX(f * (1 << _cie_interpolate_bits),\
  856.               (n) << _cie_interpolate_bits, itemp)
  857. #else
  858.         int m = pcrd->RenderTable.lookup.m;
  859.         int k =
  860.             (c == 0 ? 1 : c == 1 ?
  861.              m * pcrd->RenderTable.lookup.dims[2] : m);
  862. #  define SCALED_INDEX(f, n, itemp)\
  863.      (RESTRICTED_INDEX(f, n, itemp) * k)
  864. #endif
  865.         const gs_range *prange = pcrd->RangeABC.ranges + c;
  866.         double scale = (n - 1) / (prange->rmax - prange->rmin);
  867.  
  868.         for (i = 0; i < gx_cie_cache_size; ++i) {
  869.             float v =
  870.             (pcache->floats.values[i] - prange->rmin) * scale
  871. #ifndef CIE_RENDER_TABLE_INTERPOLATE
  872.             + 0.5
  873. #endif
  874.             ;
  875.             int itemp;
  876.  
  877.             if_debug5('c',
  878.                   "[c]cache[%d][%d] = %g => %g => %d\n",
  879.                   c, i, pcache->floats.values[i], v,
  880.                   SCALED_INDEX(v, n, itemp));
  881.             pcache->fixeds.ints.values[i] =
  882.             SCALED_INDEX(v, n, itemp);
  883.         }
  884.         pcache->fixeds.ints.params = pcache->floats.params;
  885.         pcache->fixeds.ints.params.is_identity = false;
  886. #undef SCALED_INDEX
  887.         }
  888.     }
  889.     /* Fold the scaling of the EncodeABC cache index */
  890.     /* into MatrixABC. */
  891. #define MABC(i, t)\
  892.   f = pcrd->caches.EncodeABC[i].floats.params.factor;\
  893.   pcrd->MatrixABCEncode.cu.t *= f;\
  894.   pcrd->MatrixABCEncode.cv.t *= f;\
  895.   pcrd->MatrixABCEncode.cw.t *= f;\
  896.   pcrd->EncodeABC_base[i] =\
  897.     float2cie_cached(pcrd->caches.EncodeABC[i].floats.params.base * f)
  898.     MABC(0, u);
  899.     MABC(1, v);
  900.     MABC(2, w);
  901. #undef MABC
  902.     pcrd->MatrixABCEncode.is_identity = 0;
  903.     }
  904.     cie_cache_mult3(pcrd->caches.EncodeLMN, &pcrd->MatrixABCEncode);
  905.     pcrd->status = CIE_RENDER_STATUS_COMPLETED;
  906.     return 0;
  907. }
  908.  
  909. /* Apply a range restriction to a cache. */
  910. private void
  911. cie_cache_restrict(cie_cache_floats * pcache, const gs_range * prange)
  912. {
  913.     int i;
  914.  
  915.     for (i = 0; i < gx_cie_cache_size; i++) {
  916.     float v = pcache->values[i];
  917.  
  918.     if (v < prange->rmin)
  919.         pcache->values[i] = prange->rmin;
  920.     else if (v > prange->rmax)
  921.         pcache->values[i] = prange->rmax;
  922.     }
  923. }
  924.  
  925. /* Convert a cache from floats to fracs. */
  926. /* Note that the two may be aliased. */
  927. void
  928. gs_cie_cache_to_fracs(const cie_cache_floats *pfloats, cie_cache_fracs *pfracs)
  929. {
  930.     int i;
  931.  
  932.     /* Loop from bottom to top so that we don't */
  933.     /* overwrite elements before they're used. */
  934.     for (i = 0; i < gx_cie_cache_size; ++i)
  935.     pfracs->values[i] = float2frac(pfloats->values[i]);
  936.     pfracs->params = pfloats->params;
  937. }
  938.  
  939. /* ------ Fill in the joint cache ------ */
  940.  
  941. /* If the current color space is a CIE space, or has a CIE base space, */
  942. /* return a pointer to the common part of the space; otherwise return 0. */
  943. private const gs_cie_common *
  944. cie_cs_common_abc(const gs_color_space *pcs_orig, const gs_cie_abc **ppabc)
  945. {
  946.     const gs_color_space *pcs = pcs_orig;
  947.  
  948.     *ppabc = 0;
  949.     do {
  950.         switch (pcs->type->index) {
  951.     case gs_color_space_index_CIEDEF:
  952.         *ppabc = (const gs_cie_abc *)pcs->params.def;
  953.         return &pcs->params.def->common;
  954.     case gs_color_space_index_CIEDEFG:
  955.         *ppabc = (const gs_cie_abc *)pcs->params.defg;
  956.         return &pcs->params.defg->common;
  957.     case gs_color_space_index_CIEABC:
  958.         *ppabc = pcs->params.abc;
  959.         return &pcs->params.abc->common;
  960.     case gs_color_space_index_CIEA:
  961.         return &pcs->params.a->common;
  962.     default:
  963.             pcs = gs_cspace_base_space(pcs);
  964.             break;
  965.         }
  966.     } while (pcs != 0);
  967.  
  968.     return 0;
  969. }
  970. const gs_cie_common *
  971. gs_cie_cs_common(const gs_state * pgs)
  972. {
  973.     const gs_cie_abc *ignore_pabc;
  974.  
  975.     return cie_cs_common_abc(pgs->color_space, &ignore_pabc);
  976. }
  977.  
  978. /*
  979.  * Mark the joint caches as needing completion.  This is done lazily,
  980.  * when a color is being mapped.  However, make sure the joint caches
  981.  * exist now.
  982.  */
  983. int
  984. gs_cie_cs_complete(gs_state * pgs, bool init)
  985. {
  986.     gx_cie_joint_caches *pjc = gx_currentciecaches(pgs);
  987.  
  988.     if (pjc == 0)
  989.     return_error(gs_error_VMerror);
  990.     pjc->status = (init ? CIE_JC_STATUS_BUILT : CIE_JC_STATUS_INITED);
  991.     return 0;
  992. }
  993. /* Actually complete the joint caches. */
  994. int
  995. gs_cie_jc_complete(const gs_imager_state *pis, const gs_color_space *pcs)
  996. {
  997.     const gs_cie_abc *pabc;
  998.     const gs_cie_common *common = cie_cs_common_abc(pcs, &pabc);
  999.     gs_cie_render *pcrd = pis->cie_render;
  1000.     gx_cie_joint_caches *pjc = pis->cie_joint_caches;
  1001.  
  1002.     if (pjc->cspace_id == pcs->id &&
  1003.     pjc->render_id == pcrd->id)
  1004.     pjc->status = pjc->id_status;
  1005.     switch (pjc->status) {
  1006.     case CIE_JC_STATUS_BUILT: {
  1007.     int code = cie_joint_caches_init(pjc, common, pcrd);
  1008.  
  1009.     if (code < 0)
  1010.         return code;
  1011.     }
  1012.     /* falls through */
  1013.     case CIE_JC_STATUS_INITED:
  1014.     cie_joint_caches_complete(pjc, common, pabc, pcrd);
  1015.     pjc->cspace_id = pcs->id;
  1016.     pjc->render_id = pcrd->id;
  1017.     pjc->id_status = pjc->status = CIE_JC_STATUS_COMPLETED;
  1018.     /* falls through */
  1019.     case CIE_JC_STATUS_COMPLETED:
  1020.     break;
  1021.     }
  1022.     return 0;
  1023. }
  1024.  
  1025. /*
  1026.  * Compute the source and destination WhitePoint and BlackPoint for
  1027.  * the TransformPQR procedure.
  1028.  */
  1029. int 
  1030. gs_cie_compute_points_sd(gx_cie_joint_caches *pjc,
  1031.              const gs_cie_common * pcie,
  1032.              const gs_cie_render * pcrd)
  1033. {
  1034.     gs_cie_wbsd *pwbsd = &pjc->points_sd;
  1035.  
  1036.     pwbsd->ws.xyz = pcie->points.WhitePoint;
  1037.     cie_mult3(&pwbsd->ws.xyz, &pcrd->MatrixPQR, &pwbsd->ws.pqr);
  1038.     pwbsd->bs.xyz = pcie->points.BlackPoint;
  1039.     cie_mult3(&pwbsd->bs.xyz, &pcrd->MatrixPQR, &pwbsd->bs.pqr);
  1040.     pwbsd->wd.xyz = pcrd->points.WhitePoint;
  1041.     pwbsd->wd.pqr = pcrd->wdpqr;
  1042.     pwbsd->bd.xyz = pcrd->points.BlackPoint;
  1043.     pwbsd->bd.pqr = pcrd->bdpqr;
  1044.     return 0;
  1045. }
  1046.  
  1047. /*
  1048.  * Sample the TransformPQR procedure for the joint caches.
  1049.  * This routine is idempotent.
  1050.  */
  1051. private int
  1052. cie_joint_caches_init(gx_cie_joint_caches * pjc,
  1053.               const gs_cie_common * pcie,
  1054.               gs_cie_render * pcrd)
  1055. {
  1056.     bool is_identity;
  1057.     int j;
  1058.  
  1059.     gs_cie_compute_points_sd(pjc, pcie, pcrd);
  1060.     /*
  1061.      * If a client pre-loaded the cache, we can't adjust the range.
  1062.      * ****** WRONG ******
  1063.      */
  1064.     if (pcrd->TransformPQR.proc == TransformPQR_from_cache.proc)
  1065.     return 0;
  1066.     is_identity = pcrd->TransformPQR.proc == TransformPQR_default.proc;
  1067.     for (j = 0; j < 3; j++) {
  1068.     int i;
  1069.     gs_for_loop_params lp;
  1070.  
  1071.     gs_cie_cache_init(&pjc->TransformPQR[j].floats.params, &lp,
  1072.               &pcrd->RangePQR.ranges[j], "TransformPQR");
  1073.     for (i = 0; i < gx_cie_cache_size; lp.init += lp.step, i++) {
  1074.         float out;
  1075.         int code =
  1076.         (*pcrd->TransformPQR.proc)(j, lp.init, &pjc->points_sd,
  1077.                        pcrd, &out);
  1078.  
  1079.         if (code < 0)
  1080.         return code;
  1081.         pjc->TransformPQR[j].floats.values[i] = out;
  1082.         if_debug4('C', "[C]TransformPQR[%d,%d] = %g => %g\n",
  1083.               j, i, lp.init, out);
  1084.     }
  1085.     pjc->TransformPQR[j].floats.params.is_identity = is_identity;
  1086.     }
  1087.     return 0;
  1088. }
  1089.  
  1090. /*
  1091.  * Complete the loading of the joint caches.
  1092.  * This routine is idempotent.
  1093.  */
  1094. private void
  1095. cie_joint_caches_complete(gx_cie_joint_caches * pjc,
  1096.               const gs_cie_common * pcie,
  1097.               const gs_cie_abc * pabc /* NULL if CIEA */,
  1098.               const gs_cie_render * pcrd)
  1099. {
  1100.     gs_matrix3 mat3, mat2, mat1;
  1101.     gs_matrix3 MatrixLMN_PQR;
  1102.     int j;
  1103.  
  1104.     /*
  1105.      * We number the pipeline steps as follows:
  1106.      *   1 - DecodeABC/MatrixABC
  1107.      *   2 - DecodeLMN/MatrixLMN/MatrixPQR
  1108.      *   3 - TransformPQR/MatrixPQR'/MatrixLMN
  1109.      *   4 - EncodeLMN/MatrixABC
  1110.      *   5 - EncodeABC, RenderTable (we don't do anything with this here)
  1111.      * We work from back to front, combining steps where possible.
  1112.      * Currently we only combine steps if a procedure is the identity
  1113.      * transform, but we could do it whenever the procedure is linear.
  1114.      * A project for another day....
  1115.      */
  1116.     /* Step 4 */
  1117.     if (pcrd->caches.EncodeLMN[0].floats.params.is_identity &&
  1118.     pcrd->caches.EncodeLMN[1].floats.params.is_identity &&
  1119.     pcrd->caches.EncodeLMN[2].floats.params.is_identity) {
  1120.     /* Fold step 4 into step 3. */
  1121.     cie_matrix_mult3(&pcrd->MatrixABCEncode, &pcrd->MatrixPQR_inverse_LMN,
  1122.              &mat3);
  1123.     pjc->skipEncodeLMN = true;
  1124.     } else {
  1125.     mat3 = pcrd->MatrixPQR_inverse_LMN;
  1126.     pjc->skipEncodeLMN = false;
  1127.     }
  1128.     /* Step 3 */
  1129.     cache3_set_linear(pjc->TransformPQR);
  1130.     cie_matrix_mult3(&pcrd->MatrixPQR, &pcie->MatrixLMN,
  1131.              &MatrixLMN_PQR);
  1132.     if (pjc->TransformPQR[0].floats.params.is_identity &
  1133.     pjc->TransformPQR[1].floats.params.is_identity &
  1134.     pjc->TransformPQR[2].floats.params.is_identity) {
  1135.     /* Fold step 3 into step 2. */
  1136.     cie_matrix_mult3(&mat3, &MatrixLMN_PQR, &mat2);
  1137.     pjc->skipPQR = true;
  1138.     } else {
  1139.     mat2 = MatrixLMN_PQR;
  1140.     for (j = 0; j < 3; j++) {
  1141.         cie_cache_restrict(&pjc->TransformPQR[j].floats,
  1142.                    &pcrd->RangePQR.ranges[j]);
  1143.     }
  1144.     cie_cache_mult3(pjc->TransformPQR, &mat3);
  1145.     pjc->skipPQR = false;
  1146.     }
  1147.     /* Steps 2 & 1 */
  1148.     if (pcie->caches.DecodeLMN[0].floats.params.is_identity &
  1149.     pcie->caches.DecodeLMN[1].floats.params.is_identity &
  1150.     pcie->caches.DecodeLMN[2].floats.params.is_identity) {
  1151.     if (!pabc) {
  1152.         pjc->skipDecodeLMN = mat2.is_identity;
  1153.         pjc->skipDecodeABC = false;
  1154.         if (!pjc->skipDecodeLMN)
  1155.         for (j = 0; j < 3; j++) {
  1156.             cie_cache_mult(&pjc->DecodeLMN[j], &mat2.cu + j,
  1157.                    &pcie->caches.DecodeLMN[j].floats);
  1158.         }
  1159.     } else {
  1160.         /*
  1161.          * Fold step 2 into step 1.  This is a little different because
  1162.          * the data for step 1 are in the color space structure.
  1163.          */
  1164.         cie_matrix_mult3(&mat2, &pabc->MatrixABC, &mat1);
  1165.         for (j = 0; j < 3; j++) {
  1166.         cie_cache_mult(&pjc->DecodeLMN[j], &mat1.cu + j,
  1167.                    &pabc->caches.DecodeABC[j].floats);
  1168.         }
  1169.         pjc->skipDecodeLMN = false;
  1170.         pjc->skipDecodeABC = true;
  1171.     }
  1172.     } else {
  1173.     for (j = 0; j < 3; j++) {
  1174.         cie_cache_mult(&pjc->DecodeLMN[j], &mat2.cu + j,
  1175.                &pcie->caches.DecodeLMN[j].floats);
  1176.     }
  1177.     pjc->skipDecodeLMN = false;
  1178.     pjc->skipDecodeABC = pabc != 0 && pabc->caches.skipABC;
  1179.     }
  1180. }
  1181.  
  1182. /* ================ Utilities ================ */
  1183.  
  1184. /* Multiply a vector by a matrix. */
  1185. /* Note that we are computing M * V where v is a column vector. */
  1186. private void
  1187. cie_mult3(const gs_vector3 * in, register const gs_matrix3 * mat,
  1188.       gs_vector3 * out)
  1189. {
  1190.     if_debug_vector3("[c]mult", in);
  1191.     if_debug_matrix3("    *", mat);
  1192.     {
  1193.     float u = in->u, v = in->v, w = in->w;
  1194.  
  1195.     out->u = (u * mat->cu.u) + (v * mat->cv.u) + (w * mat->cw.u);
  1196.     out->v = (u * mat->cu.v) + (v * mat->cv.v) + (w * mat->cw.v);
  1197.     out->w = (u * mat->cu.w) + (v * mat->cv.w) + (w * mat->cw.w);
  1198.     }
  1199.     if_debug_vector3("    =", out);
  1200. }
  1201.  
  1202. /*
  1203.  * Multiply two matrices.  Note that the composition of the transformations
  1204.  * M1 followed by M2 is M2 * M1, not M1 * M2.  (See gscie.h for details.)
  1205.  */
  1206. private void
  1207. cie_matrix_mult3(const gs_matrix3 *ma, const gs_matrix3 *mb, gs_matrix3 *mc)
  1208. {
  1209.     gs_matrix3 mprod;
  1210.     gs_matrix3 *mp = (mc == ma || mc == mb ? &mprod : mc);
  1211.  
  1212.     if_debug_matrix3("[c]matrix_mult", ma);
  1213.     if_debug_matrix3("             *", mb);
  1214.     cie_mult3(&mb->cu, ma, &mp->cu);
  1215.     cie_mult3(&mb->cv, ma, &mp->cv);
  1216.     cie_mult3(&mb->cw, ma, &mp->cw);
  1217.     cie_matrix_init(mp);
  1218.     if_debug_matrix3("             =", mp);
  1219.     if (mp != mc)
  1220.     *mc = *mp;
  1221. }
  1222.  
  1223. /* Invert a matrix. */
  1224. /* The output must not be an alias for the input. */
  1225. private void
  1226. cie_invert3(const gs_matrix3 *in, gs_matrix3 *out)
  1227. {    /* This is a brute force algorithm; maybe there are better. */
  1228.     /* We label the array elements */
  1229.     /*   [ A B C ]   */
  1230.     /*   [ D E F ]   */
  1231.     /*   [ G H I ]   */
  1232. #define A cu.u
  1233. #define B cv.u
  1234. #define C cw.u
  1235. #define D cu.v
  1236. #define E cv.v
  1237. #define F cw.v
  1238. #define G cu.w
  1239. #define H cv.w
  1240. #define I cw.w
  1241.     double coA = in->E * in->I - in->F * in->H;
  1242.     double coB = in->F * in->G - in->D * in->I;
  1243.     double coC = in->D * in->H - in->E * in->G;
  1244.     double det = in->A * coA + in->B * coB + in->C * coC;
  1245.  
  1246.     if_debug_matrix3("[c]invert", in);
  1247.     out->A = coA / det;
  1248.     out->D = coB / det;
  1249.     out->G = coC / det;
  1250.     out->B = (in->C * in->H - in->B * in->I) / det;
  1251.     out->E = (in->A * in->I - in->C * in->G) / det;
  1252.     out->H = (in->B * in->G - in->A * in->H) / det;
  1253.     out->C = (in->B * in->F - in->C * in->E) / det;
  1254.     out->F = (in->C * in->D - in->A * in->F) / det;
  1255.     out->I = (in->A * in->E - in->B * in->D) / det;
  1256.     if_debug_matrix3("        =", out);
  1257. #undef A
  1258. #undef B
  1259. #undef C
  1260. #undef D
  1261. #undef E
  1262. #undef F
  1263. #undef G
  1264. #undef H
  1265. #undef I
  1266.     out->is_identity = in->is_identity;
  1267. }
  1268.  
  1269. /* Set the is_identity flag that accelerates multiplication. */
  1270. private void
  1271. cie_matrix_init(register gs_matrix3 * mat)
  1272. {
  1273.     mat->is_identity =
  1274.     mat->cu.u == 1.0 && is_fzero2(mat->cu.v, mat->cu.w) &&
  1275.     mat->cv.v == 1.0 && is_fzero2(mat->cv.u, mat->cv.w) &&
  1276.     mat->cw.w == 1.0 && is_fzero2(mat->cw.u, mat->cw.v);
  1277. }
  1278.